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A FEASIBILITY STUDY 
. of: A 

3-D FINITE ELEMENT SOLUTION SCHEME 

FOR 

AEROENGINE DUCT ACOUSTICS 
By 


A. L. Abrahamson 
Wyle Laboratories 

INTRODUCTION 


Development of two-dimensional (2-D) finite element solution schemes for 
aeroengine duct acoustics was initiated in the mid-1970s and has achieved con- 
siderable success. Analysis of acoustic propagation in axisymmetric ducts with 
nonuniform flow was made possible^ ^ and schemes for increasing computational 

efficiency have progressed to the stage where repeated analyses for aero- 

( 2 ) 

engine duct liner optimization are possible. 


The satisfactory results with 2-D models and the ability of segmented 
duct liners to give good sound attenuation characteristics indicated a potential 
advantage from the development of three-dimensional (3-D) finite element models. 


The expected advantage was based on the rationale that since axially 
segmented acoustic duct liners have been demonstrated to provide improved 
sound attenuation characteristics over uniform liners, there is a significant 
likelihood that a combination of axially and circumferentially segmented liners 
will provide still better sound attenuation characteristics. In order to model 
the effect of this type of liner segmentation, it is necessary to discard the 
axisymmetric assumption and develop a fully 3-D model. 

With currently available computers, this was anticipated to be ah 
ambitious task. However, the use of different element formulation techniques 
and different linear equation solution algorithms to those used in the 2-D 
analyses indicated that the task might be possible. 



SYMBOLS 


A 

b 

B 

c 

I 

L,M,N 

N 

0 

N. 

1 

P 

u,v,w 


u,v,w 

X 

x,y,z 

n 

to 

5,n> s 
p 


Global matrix 

Rigbt-band side of global matrix equation 
Matrix bandwidth 
Speed of sound 
Identity of matrix 

Number of finite elements along coordinate directions 
of rectangular mesh 

Global matrix order 

Finite element functional 

Acoustic pressure 

Acoustic velocities along coordinate directions 

Mean aerodynamic velocities along coordinate directions 

Solution vector 

Cartesian coordinates 

Over-relaxation parameter 

Acoustic frequency 

Nodal parameter 

Element local coordinates 

Mean density 


CHARACTERISTICS OF THE PROBLEM 
Order 

Numerical solution of the differential equations describing sound pro- 
pagation in an aeroengine inlet requires discretization of the space within 
which the solution is required. The order of the discretization is dependent 
upon the spatial rate of change of the variables, which in turn is dependent 
upon acoustic source frequency and geometric and aerodynamic parameters. 
Precise calculation of the required order of discretization is not possible since 
the problem solution is not known a priori. Approximate calculations of 
required linear finite element discretization meshes for 2-D and 3-D analyses 
of typical aeroengine inlets, based upon a linear one-dimensional analytical 
model, are shown in Table I. 
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Table I 


Approximate Discretization Order for Solution of 
Sound Propagation Equations in an Aeroengine Inlet 

Diameter 

(Meters) 

Length 

(Meters) 

Mach 

No. 

Frequency 

(kHz) 


2-D 

Solution 


3-D Solution || 

M 

N 

Matrix 

Order 

(N o ) 

Matrix 

Band- 

Width 

(B) 

1 

M 

N 

Matrix 

Order 

(N q ) 

Matrix 
Band- 
Width 
(B) 1 

1.5 

1.5 

0 

1 

20 

40 

6,400 

320 

40 

40 

40 

512,000 

51,200 




2 

40 

80 

25,600 

640 

80 

80 

80 

4 x 10 6 

2 x 10 5 




3 

60 

120 

57,600 

960 

120 

120 

120 

13.8 x 10 6 

4.6 x 10 5 



-.3 

1 

20 

60 

9,600 

320 

40 

40 

60 

768,000 

51,200 




2 

40 

120 

38,400 

640 

80 

80 

120 

6.1 x 10 6 

2 x 10 5 




3 

60 

180 

86,400 

960 

120 

120 

180 

20.7 x 10 6 

4.6 x 10 5 

2.2 

2.2 

0 

1 

30 

. 60 

14,400 

480 

60 

60 

60 

1.7 x 10 6 

1.15 x 10 5 




2 

60 

120 

57,600 

960 

120 

120 

120 

13.8 x 10 6 

4.6 x 10 5 




3 

90 

180 

129,600 

1,440 

180 

180 

180 

46.6 x 18 

1.04 x 10 6 



-.3 

1 

30 

90 

21,600 

480 

60 

60 

90 

2.6 x 10 6 

1.15 x 10 5 




2 

60 

180 

86,400 

960 

120 

120 

180 

20.7 x 10 6 

4.6 x 10 5 




3 

90 

270 

194,400 

1,440 

180 

180 

270 

70 x 10 6 

1.04 x 10 6 


Note: Assumes linear elements, approximately 10 elements /wavelength, 8 variables /node. 

2D Solution - N = 8MN , B = 2(8M). 3D Solution - N = 8LMN , B = 4(8LM). 
o o 


The global matrix orders (N 0 ) and bandwidths (B) resulting from these 
discretization meshes are also shown in Table I. For the 2-D model, matrix 
orders range from 6,400 to 194,400; while for the 3-D model, orders range 
from 0.5 million to 70 million. It should be emphasized that the wide range 
spanned by these numerical values represent relatively small changes in duct 
physical dimensions, mean-flow characteristics and source frequency. Matrix 
bandwidths span a similarly wide range, from 320 to 1,440 in the 2-D case, 
and 51,200 to 1.04 million in the 3-D case. 


Computational Effort 

Consider matrices of the dimensions described above relative to the 
computational effort required to generate and solve them. Absolute eval- 
uation of computational effort is extremely complicated and is dependent upon 
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computer and compiler characteristics as well as on programming techniques. 

The approach taken in this determination is based upon an approximate eval- 
uation of operation counts followed by scaling of central processor unit (CPU) 
time taken from a baseline optimized FORTRAN program^ ^ on a CDC Cyber 175 
computer. 

The basic finite element process involves two distinct phases, i.e., 
assembling the global matrix, and solving the global matrix. 

In the case of boundary condition optimization, this process is augmented 
by a third phase; specifically, use of the current solution to generate an 
increment in boundary conditions. The complete process is repeated until 
convergence to an optimum set of boundary conditions is obtained. 

By defining an operation as one addition followed by one multiplication, 
we may estimate operation counts as follows: 

• Global matrix assembly - For each element, these procedures 
involve numerical integration of the finite element approximation to 
the differential equations, over the domain. The procedures are 
proportional to the number of finite elements and are thus pro- 
portional to the number of degrees of freedom (N 0 ) . The constant 
of proportionality is typically of the order 10 2 for 2-D models and 
10 3 for 3-D models. 

• Global matrix solution - Solution of a full matrix is most commonly 
achieved by a direct method such as Gaussian elimination which 
requires approximately 1/3 N 3 operations. From Table I, however, 
it may be seen that the global matrix is not full, but possesses a 
bandwidth at least one order of magnitude (sometimes 2 orders of 
magnitude in the 2-D case) less than the glob al matrix order. In 
these circumstances, Gaussian elimination may be carried out in 
approximately i N q B 2 op 61 "^ 0113 * 
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• Boundary condition increment generation - Any of the standard 

optimization methods may be used; Davidon-Fletcher-Powell, Simplex, 
or Broyden-Fletcher-Shanno, for example. The computational effort 
in implementing any of these methods is negligible and the impact 
of this phase of the process may be ignored. It is important to 
note, however, that the efficiency of the method used has a vital 
bearing on overall computational time. The most efficient method 
requires the minimum iterations for convergence to an optimum liner 
configuration . 

The operation counts given above may now be scaled by CPU times 
(Table II) from a baseline 2-D finite element model of sound propagation in 
an aeroengine duct. This model uses a numerically integrated bilinear 
isoparametric serendipity element and a block tridiagonal solution scheme 
specifically coded for problem sparsity. Figures 1 and 2 show that a linear 
relationship does hold between the theoretical operation counts given above 
the actual CPU times for assembly and solution phases. 


Table II 


Computing Times for Baseline Model 



Mat 



CPU Times for 

1/4 N B 2 
o 

(Operations) 

M 

N 

Order 

(N 0 > 


B 

Assembly 

(Sec) 


40 

110 

36,408 


656 

330 

2,060 . 

. 39 x 10 10 

25 

59 

12,480 


416 

98 

299 

. 55 x 10 9 

12 

30 

3,224 


208 

25 

35 

.35 x 10 8 

20 

200 

33,768 


336 

304 

603 

. 95 x 10 9 


Estimated computer times for the aeroengine inlet examples given in 
Table I, derived from the relations: 
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CPU I . , = 

| assembly 

CPU |solution = 
are given in Table III. 


(.94 x l(f 2 ) N 0 
(1.33 x 10~ 7 )N o B 


(secs) 

2 (secs) 


Table III 


Estimated CPU Times for Solution of Sound Propagation 
Equations in Typical Aeroengine Inlets 




CPU 

Times 


Diameter 

(Meters) 

Length 

(Meters) 

Mach 

No. 

Frequency 

(kHz) 

2- 

.D 

3. 

.D 

Assembly 

(Sec) 

Solution 

(Min) 

Assembly 

(Min) 

msmM 

1.5 

1.5 

0 

1 

6 

1.5 

8 

4.9 x 10 4 




2 

24 

23 

62 

5.9 x 10 6 




3 

54 

118 

216 




-.3 

1 

9 

2.2 

12 





2 

36 

35 

95 





3 

81 

176 

324 


2.2 

2.2 

0 

1 

14 

7.4 

27 





2 

54 

118 

216 





3 

122 

596 

730 




-.3 

1 

20 

11 

41 





2 

81 

177 

324 

/ 




3 

183 

894 

1,097 



NOTE : Assumes linear finite elements, approximately 10 elements/wavelength, 

8 variables /node. 
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Figure 1. - Plot of theoretical operation counts (N ) versus 

o 

actual matrix assembly times for baseline model. 
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MATRIX SOLUTION 



C. 

Figure 2. - Plot of theoretical operation counts (N^B ) versus 
actual matrix solution times for baseline model. 



8 





CPU time estimates from Table III show that global matrix assembly- 
times are substantially less than solution times. For a 2-D analysis, assembly 
times are measured in seconds while solution times range from a few minutes 
to several hours. From a cursory examination of solution times, it is 
evidently pointless to ever contemplate a 3-D analysis using the techniques 
of the baseline model since CPU times begin at 50,000 hours. 

A 2-D analysis of a single case spans the range from 1.5 minutes, for 
a 1-kHz source and a zero Mach number in a duct 1. 5 m in diameter and 1.5 m 
long, to 15 hours for a 3-kHz source with a flow Mach number of -0.3 in a duct 
2.2 m in diameter and 2.2 m long. 

Suppose we take one hour as a reasonable upper limit of CPU time per 
case, then a 2-D analysis of the 1.5-m duct,' up to a source frequency between 
2 kHz and 3 kHz is feasible. Alternately, for the 2.2-m duct, the feasible 
upper frequency range lies somewhere between 1 kHz and 2 kHz. 

Consider the process of duct liner optimization. Depending upon the 
number of degrees of freedom in the duct liner specification, numerical 
optimization requires repetition of the analysis typically between 20 and 200 
times. As shown in reference 2, it is not necessary to multiply CPU times 
given in Table III by this factor, since alternate direct decomposition techniques 
such as the boundary exclusive decomposition method are available. 


Alternate Approaches for 3-D Analysis 

Since the methods used in 2-D models are clearly inadequate for a 
3-D analysis, new approaches are necessary. It has been demonstrated that 
the vast increases in computer time from a 2-D to a 3-D analysis confe from 
attempting a direct solution to a matrix equation whose bandwidth has 
increased by two or three orders of magnitude. 

The alternative to a direct method is an iterative method. Iterative 
solution algorithms such as Jacobi, Gauss-Seidell and Successive Over- 
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Relaxation (SOR) differ from Gaussian elimination in. one important respect. 
Whereas, Gaussian elimination generates "fill," i.e., large quantitites of 
intermediate numbers which all require processing, these algorithms always 
work only on the original matrix. 

To provide an indication of the potential savings which may be gained 
from implementation of these algorithms, consider the original number of 
nonzeroes per row in the global matrix for the examples given in Table I. 

For all 2-D discretizations, these number just 72; while for all 3-D discre- 
tizations, 216 nonzeroes exist per row. Comparing these numbers with the 
matrix bandwidths shown in Table I, the original degree of sparsity is evident. 
In the process of Gaussian elimination, the entire bandwidth is filled in the 
solution process. 

Two problems exist with iterative solutions, however. They are as 
follows : 

• A positive definite global matrix is required for convergence to be 
guaranteed. 

• The number of iterations required for adequate convergence may 
be very small or very large depending upon the characteristics 
of the matrix. 

The Galerkin method of formulating finite element equations does not give a 
positive definite global matrix. The least squares method, however, does 
give a symmetric positive definite global matrix but frequently the least 
squares process yields equations that are poorly conditioned. 

i 

Estimated CPU times for one iteration of each of the typical inlet examples 
are given in Table IV. 
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Table IV 


Estimated CPU Times /Iteration for Iterative Solution 
of Sound Propagation equations in an' Aeroengine Inlet 


Diameter 

(Meters) 

Length 

(Meters) 

Mach 

No. 

Frequency 

(kHz) 

2-D 
(Min. ) 

3-D 

(Min) 

1.5 

' 1.5 

0 

1 

1.02 x 10~ 3 

.25 




2 

3.99 x 10" 3 

1.9 




3 

9.09 x 10" 3 

6.6 



-.3 

1 

1.53 x 10" 3 

.37 




2 

6.12. x 10~ 3 

2.9 




3 

1.37 x 10“ 2 

9.9 

2.2 

2.2 

0 

1 

2.30 x 10~ 3 

.81 




• 

2 

9.19 x 10~ 3 

6. 6 




3 

2.07 x 10 -2 

22 



-•3 

1 

3.44 x 10~ 3 

1.24 




2 

1.38 x 10~ 2 

9.9 




3 

3.10 x 10 _2 

33.5 


The critical question regarding the feasibility of a 3-D analysis hinges on 
the convergence rate. That is, how many iterations are required for con- 
vergence. Since convergence rate is problem dependent, the only way to 
determine it is by first assembling the global matrix equation. Thus., for 
practical purposes, the only way to determine the feasibility of the approach 
is by trying it. 


DERIVATION 


Derivation of the equations governing linearized acoustic motion for the 
2-D case, is shown in reference 1. Extension of this derivation to 3-D is 
readily performed and will not be given here. The four equations are: 
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, au,-9u,-3u,-au, au , 3u 

10 U + U — + U -r + V - + W - — + V -r— + w— 

9x 9 x 9y 9 z 9y 3z 


, p /- 3u , - 3u , - 3u N , 1 9p _ n 
jf* ^ sx" v 3y + ' > 'aT jj” 3x ' 0 


(l) 


9v , 9v , - 9 v , - 3 v , - 9v , 9v 

u p— + iov + v — + u tt — + v -r — + w t— + w — 

3x 9y 9x 3y 3z 3z 


, p /- 9v , - 3v - 9v N , 1 9p n 

+ :t 2 (ux— + Vt— + w— ) + — = 0 

pc 2 9x 3y 9z - 9y 


( 2 ) 


9w , 9w . , 3w , - 9w , - 3w , - 9w 

u __ + v— - +01W + Wr— +Ur— + Vt— + Wr— 

3x 3y 3z 9x 3y 3z 


P / - 3w . - 9w , - 9w N , 1 9p n 
+ -*- (u — + v— + w— -) + — = 0 

o „ 2 9x 9y 3z - 9z 
h<- D 


( 3 ) 


yu_ i£ + \ + ( 1 iP-f 5! i + / I Ip . 3w x 

- 9x 9x ' - 9y 9y ' ' - dz 9z ' 

P P 3 3 P 


+ — [up + u^ + v^ + w — + p(— 


pc' 


gx 9y 9Z 9X 


+ 


9 v 

9y 


+ 


3w 

9z 


) ] = 0 


( 4 ) 


where (u,v,w) and (u,v,w) are the (x,y,z) components of acoustic and 
mean-flow velocities, respectively, p is the acoustic pressure, and p the mean 
density. Acoustic source frequency is u and c is the local speed of sound. 
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(3) 

A 3-D cubic serendiptity element was selected for the analysis, This 
element shown in Figure 3 has a total of 32 nodes and when applied to 
solution of the above equations with four variables per node (u,v,w,p) yields 
a local element matrix of order 256. 

The principal reason for adoption of this element over a linear element was 
to take full advantage of the computational speed increases offered by a vector 
computer. With iterative solution algorithms on a vector computer (such as 
the CDC STAR 100) , vector lengths for global matrix solution using this 
element lie between 256 and 1080. These vectors are sufficiently long to 
minimize the effect of "pipeline" processor startup time on this machine. 

Consider the volume of space enclosed by the cube shown in Figure 3. 
Within this cube, a parameter <|>(z,y,z) is completely defined by its values 
{ } e at the nodes through the relation: 

* = (N) {<f >} 6 ( Ji ) 

= (i V n 2 ’ * • * N 3 2 ) j : 2 | 

* * 32 * ( 5 ) 

Functionals X N. where i ■= 1, 4, 9, 12, 21, 24, 29, 32 take the form: 
lN i = U (1 + (1 + n o ) (1 + 5 0 > {9U 2 +C 2 )-19} 

Functionals 2 N. where i — 2, 3, 10, 11, 22, 23, 30, 31 take the form: 

2 N = _9 

i 64 (1 -t, 1 ) (1 + % „) (1 + O 0 ) (1 +C 0 ) 

f 

/ 

Functionals where i — 5, 6, 7, 8, 25, 26, 27, 28 take the form: 

3N i = U (1 ■ n2) U + 9 °0> O * s 0 ) <1 + c„> 
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Functionals . l *N where i— 13, 14, 15, 16, 17, 18, 19, 20 take the form 
4 N. = (1 - c 2 ) n + 9? 0 ) (1+ S 0 ) (1 +. no) 


In the above expressions, 

5 0 ~ ^i’ 11 o - 11 .o'V ^ ~ ^ ^i 

where (£., n>, ?-) represent the nondimensional coordinates of the ith node 
and may assume values ±1 or ± 173. 

Suppose (u,v,w,p) and (u,v,w,p) are given by relations of the form in 
equation (5), then these functionals can be substituted into equations (1) 
through (4) . 

After carrying out an isoparametric mapping of the element equations (1) 
through (4), they are squared and summed. The variation with respect to 
each nodal variable is taken to yield a set of 256 equations. These are 
integrated numerically over the volume to give the local element matrix. A 
boundary condition matrix is compiled by squaring and summing the boundary 
condition equations integrated over element boundary surfaces and taking 
variations with respect to nodal variables. The boundary condition matrix is 
multiplied by a constant, determined by best numerical conditioning, and 
summed with the element matrix. (The Least Squares Method is described in 
detail in references 4 and 5.) 

PROGRAM ORGANIZATION 

/ 

The computer program to perform a 3-D duct acoustic analysis was 
designed in modular form, The principal modules and their functions are given 
in the Appendix. 
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SOLUTION TECHNIQUE EVALUATION 


Checkout with Direct Equation Solver 

To verify the element formulation, boundary condition insertion, and 
program coding, test cases were compiled. The first test case consisted simply 
of a single element, giving a system of equations of order 256, with straight 
sides, i.e., a cube. The x-y plane at z = 0 was specified as a source plane, 
and at a z = 1 as a radiation plane with an impedance of pc. All other walls 
were hard (impedance = °°). The source frequency was chosen to be 
sufficiently low so that only the plane-wave mode was cut on, mean flow was 
zero. 

A direct-solution library subroutine was used to solve the resulting system 
of equations. This solution coincided closely with the analytically predicted, 
plane-wave propagation. Conditioning of the set of equations was evaluated 
by investigating the range of eigenvalues. This range (a factor of 10 9 ) 
indicated poor conditioning. 'Modification of equation weighting factors in the 
element least squares formulation and the weighting factors in the boundary 
condition integrals gave only marginal improvement in conditioning . 


Jacobi and SOR Iterative Methods 

Writing the matrix equation as AX = b, we may split the global matrix 
(A) into a lower triangle (L) , a diagonal (D) and an upper triangle (U), 

A = D (L + I * U) 

/ 

Jacobi's iterative method is then written as. 


X (k+i) - (L + u) x (k) + D -l b 

where X is the previous estimate for the solution vector and 
is the new estimate. 
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Similarly, the SOR method is 

X (k+1) = (I + fiL)" 1 [{(1 -ftU} X (k) * D _1 b] 
where ft, the over-relaxation parameter, lies between 1 and 2. 

For these methods to converge, it is necessary that all eigenvalues of 
their iteration matrices, viz, 

-(L + U) and (I + ftL)" 1 {(1 - ft) I - ftU} 

( M 

respectively, should be less than unity 

Programs were written to construct these iteration matrices , obtain their 
eigenvalues, and implement the two solution schemes for the test problem 
already analyzed using the direct solver. 

For the Jacobi method, the largest eigenvalues were greater than unity 
for the test problem and as expected, solutions diverged. 

For the SOR method, the largest eigenvalues were less than, but 
extremely close to one ('0.99998) indicating slow convergence. In practice, 
convergence for this method was found to be so slow that its use was 
impractical. No noticeable improvement was found by varying ft. 


Conjugate Gradient Method 

The conjugate gradient method is unique among iterative methods since, 
in the limit, it becomes a direct method. That is, for an nth order system if 
rounding errors are ignored, the nth iteration will give the exact solution. 
Often, however, a good approximation to the exact solution is obtained for 
an iteration well below n. 

(7) 

The method was first published in 1952 but did not receive much 
attention until recently when it was shown that various preconditionings can 


17 



significantly improve its rate of convergence. In most cases, it is now 
significantly faster than SOR. 

The theory of the method is given in references 7 through 9, and only 
the mechanics of the method will be described here. 


Given a system AX = B with a positive definite nxn matrix A and a start- 
ing vector X Q compute 

r = B-AX and set p„ = r n . 
o o o 


For i = 0,1,2, - - - - n-1 compute 


a. 

l 


2 



t a 

Pi Ap. 


X. 


i+1 


= X. + 

l 


a.p. 

1*1 


r i+l 


r. 

i 


- a.Ap. 

i 


b. = 

l 



p. . . = r. , . + b.p. 

*1+1 l+l 1*1 

/ 

X^ will be the solution of the linear system if rounding errors are ne- 
glected. Often a good approximation to the solution is obtained by some X. 
for i <<n. 
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Results from application of the conjugate gradient method to. the test 
problem showed high initial promise. Without preconditioning, solution of the 
test problem (order 256) was obtained in 140 iterations. With symmetric scaling, 
the number of iterations for solution was decreased to 50, using the same 
arbitrary starting vector. 

General convergence trends were then investigated with the result that 
the total number of iterations for convergence was found to be dependent upon 
several factors : 

e The variance between the first estimate and the actual solution 
• The number of finite elements per wavelength 
© The conditioning of the global matrix 
« The order of the global matrix 

In these investigations, the Euclidean norm of the residual vector (||r^|[ 2 ) 
was used to assess convergence. For a typical first estimate, | [r | | 2 was 
of the order of 10 3 ; while for a good first estimate, | |r | | was about unity. 

02 _ 4 

Convergence appeared to reach a plateau at values of I l r j II 2 between 10 
and 10 with no further improvement obtainable regardless of the number of 
subsequent iterations. 

Dependence of the number of iterations for convergence on the accuracy 
of the first estimate was expected and is a normal characteristic of iterative 
methods. Dependence on the number of finite elements per wavelength was 
surprising. After further investigation, it was found that varying the number 
of finite elements per wavelength affects the conditioning of the global matrix, 
and that the fewer finite elements allowed per wavelength, the worse the 
conditioning became. 

It had been hoped that a firm trend could be established giving the 

number of iterations required for convergence versus the order of the global 

matrix, but it was discovered that under some circumstances the Euclidean 

norm of the residual vector ( | |r. | | ) would decrease in the usual inverse 

1 2 _ 6 
exponential manner until a low value of the error was reached (about 10 ) . 

However, the solution differed from a direct solution by unacceptably 
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large margins. Apparently the equations may become sufficiently ill- 
conditioned that solution accuracy is degraded by a significant degree. 

This ill-conditioning may be seen from the plots in Figure 4. The 
analytical plane-wave solution is compared with finite element solutions for 
one through four finite elements per acoustic wavelength. In any one con- 
figuration, the global finite element matrix equation is solved in two 
different ways: first, using a direct method and second, using the pre- 

conditioned conjugate gradient method. In the case of one finite element per 
wavelength, the two solutions not only differ substantially from the analytical 
plane-wave solution, but also differ substantially from each other. In the 
cases of two and three elements per wavelength, the finite element solutions 
fall much closer to the analytical plane-wave solution; while for the case of 
four elements per wavelength, the three solutions are coincident. These 
results indicate that the finite element global matrix equation becomes increas- 
ingly well conditioned, the more elements allowed per wavelength. 

The behavior of the norm of the residual vector | |r^ | | 2 , versus 
iteration number for the cases in Figure 4, is illustrated in Figure 5. Here, 
the error norm decreases in a similar manner to similar final values in each 
case. This may be erroneously interpreted as indicating that the conjugate 
gradient solution algorithm has solved all four problems to the same degree of 
accuracy. 

The first estimate solution vector in each case is the analytical plane- 
wave solution. Note that the error norm | |r^ | | decreases very rapidly in all 
cases and that the total number of iterations is substantially less than the 
number of degrees of freedom. 

/ 

An additional problem was the unexpectedly large number of iterations 
required to reach an acceptably low value of the error norm when the first 
estimate solution vector differed significantly from the true solution. For 
example, the case presented in Figures 6 and 7 uses a large number of 
elements per wavelength (18) which should give a well conditioned matrix. 
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Figure 6. - Comparison of finite element conjugate gradient/ 
solution with analytical plane wave solution for 
poor first estimate solution vector. 





The first estimate of the solution vector was chosen to be zero except for the 
real pressure components which were all assumed to be unity* 

The total system of equations numbered 1696. After 600 iterations the 

— 2 

norm ll r 5 || 2 was a * a va l ue °f about 10 (Figure 7) and the solution 
differed substantially from the correct solution. Only after 800 - 1000 

- 4 

iterations, the error norm had decreased to sufficiently low values of 10 
_ 5 

to 10 . At these values, the solution became an accurate copy of the 

analytical plane-wave solution (Figures 6 and 7) . 

As a final exercise, some larger finite element meshes were run in order 
to assess typical computer resources required, 
e . g . , 

Mesh - 4 x 4 x 10 elements 
Degrees of Freedom - 13,240 

Computer times on the C.D.C. Star 100 were as follows: 

Global Matrix Assembly - 312 wall clock mins. 

76 CPU clock mins 


Conjugate Gradient Solution 

- 35 wall clock secs./iteration 
5 CPU clock secs 7 iteration 

Using four elements per wavelength with a good first estimate solution 
vector, [ |r. | | decreased from 4.5 x 10 to 2.5 x 10 4 in 50 iterations. 

Although the performance of the model in this example was adequate, it 
was based on the premise of an extremely good first estimate solution vector. 
With a poor first estimate solution vector, several thousand iterations would be 
required to achieve a satisfactory solution. 

The magnitude of computer time to perform several thousand iterations 
on this problem is several hours of CPU time and several days of wall clock 
time. When consideration is also made of the fact that 4 x 4 x 10 elements is 
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not a sufficiently large mesh to analyze even fairly trivial duct acoustic pro- 
blems, continuation of the effort seems unrewarding. 

It seems probable from the problems described earlier regarding poor 
conditioning of the global matrix when too few elements were allowed per wave- 
length, that severe changes in boundary conditions, element geometry or aero- 
dynamic mean variables over a short distance, would contribute further ill- 
conditioning . 


CONCLUSIONS 

This report describes an attempt to extend finite element modeling of 
duct acoustics to fully three-dimensional problems. The difficulties which arise 
when the successful 2-D techniques of Galerkin weighted residual formulation 
with a direct matrix solver are extended to 3-D, are described. The reasons 
for believing that an iterative matrix solver might be substantially quicker than 
a direct solver are indicated. Since all known iterative solvers require 
positive definite matrices, a least squares formulation was used instead of a 
Galerkin formulation . 

The root cause of the difficulties encountered in application of iterative 
matrix solvers lies in the tendency of the least squares formulation to yield 
poorly conditioned matrices. Due to this fact, the Jacobi method diverged, 
Gauss-Seidell and SOR failed to converge at a measurable rate, while the 
conjugate gradient and even direct methods were uneven in their performance. 

Of the iterative methods tried, there is no doubt that the preconditioned 

conjugate gradient method was superior. If a formulation to yield well con- 

/ 

ditioned positive-definite matrices could be found, it seems that this solver 
would give solutions within an acceptable amount of computing effort. 

Although the goal of this work was not achieved, the steps taken and 
methods used represent a logical approach to the problem which should benefit 
future workers in this area. 
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APPENDIX 


Principal Modules and Functions 
of 

Computer Program 
AFE3D1 
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Within this structure, moduLe ASEMBLE is the most significant high-level 
module- It has the following structure. 



(Sheet 1 of 2) 


30 











Call PACK 


Module packs element matrix into 
global matrix. 



END 
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